#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Quality Evaluation                                                 #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 08/27/2012                                             #
# Last Modification: 08/27/2012                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 95
#
# MANUAL: This script reports the data processing statistics.
#
# PUBLICATION: 3D reconstruction of two-dimensional crystals: <A HREF="http://www.ncbi.nlm.nih.gov/pubmed/26093179">Arch Biochem Biophys 581, 68-77 (2015)</A>
# PUBLICATION: 3D Reconstruction from 2D crystal image and diffraction data: <A HREF="http://dx.doi.org/10.1016/S0076-6879(10)82004-X">Methods Enzymol. 482, Chapter 4, 101-129 (2010)</A>
# PUBLICATION: 2dx - Automated 3D structure reconstruction from 2D crystal data: <A HREF="http://journals.cambridge.org/action/displayAbstract?aid=1943200">Microscopy and Microanalysis 14(Suppl. 2), 1290-1291 (2008)</A>
# PUBLICATION: 2dx_merge - Data management and merging for 2D crystal images: <A HREF="http://dx.doi.org/10.1016/j.jsb.2007.09.011">J. Struct. Biol. 160(3), 375-384 (2007)</A>
#
# 
# DISPLAY: MergeResolution
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: set merge_res_limit
# DISPLAY: tempkeep
# DISPLAY: realcell
# DISPLAY: realang
# DISPLAY: ALAT
# DISPLAY: zstarrange
# DISPLAY: SYM
# DISPLAY: avrgamphsNUMBER
# DISPLAY: avrgamphsRESOL
# DISPLAY: merge_data_type
# DISPLAY: ILIST
# DISPLAY: MergeHKMAX
# DISPLAY: MergeIQMAX
# DISPLAY: merge_alsoevenodd
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set create_merged_dataset = ""
set MergeResolution = ""
set zstarwin = ""
set RESMIN = ""
set RESMAX = ""
set merge_res_limit = ""
set tempkeep = ""
set realcell = ""
set realang = ""
set ALAT = ""
set zstarrange = ""
set MergeStepSize = ""
set IBOXPHS = ""
set SYM = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set merge_data_type = ""
set ILIST = ""
set MergeHKMAX = ""
set MergeIQMAX = ""
set refbeamtilt = ""
set reftiltgeo = ""
set merge_reference = ""
set merge_ref_num = ""
set max_amp_correction = ""
set JREFL = ""
set JREFL_2D = ""
set Reflections_Unique = ""
set num_amplitudes_observed = ""
set num_phases_observed = ""
set num_reflections_fitted = ""
set num_reflections_FOM1 = ""
set num_reflections_FOM50 = ""
set overall_R_factor = ""
set overall_phase_residual = ""
set overall_weighted_R_factor = ""
set overall_weighted_phase_residual = ""
set overall_phase_residual_2D = ""
set overall_weighted_phase_residual_2D = ""
set zstarrange_real = ""
set latline_algo = ""
set merge_alsoevenodd = ""
#
#$end_vars
#
setenv OMP_NUM_THREADS 8
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set date = `date`
echo date = ${date}
#
set scriptname = 2dx_QualityEvaluation
#
set merge_modus = "3D"
#
\rm -f LOGS/${scriptname}.results
#
# The following is to make sure that for the next "Import Images", the default is correctly initialized.
set initialization_reset = "y"
set initialization_executable = "y"
echo "set initialization_reset = ${initialization_reset}" >> LOGS/${scriptname}.results
echo "set initialization_executable = ${initialization_executable}" >> LOGS/${scriptname}.results
#
echo "<<@progress: 1>>"
#
#
#############################################################################
${proc_2dx}/linblock "Sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
############################################################################# 
${proc_2dx}/lin "2dx_merge_makedirs - to create all required subdirectories"
#############################################################################
#
source ${proc_2dx}/2dx_merge_makedirs
#
echo "<<@progress: 5>"
#
############################################################################# 
${proc_2dx}/lin "Gathering statistics"
#############################################################################
#
set tmp = `wc -l 2dx_merge_dirfile.dat`
set tmp_num_lines = `echo ${tmp} | awk '{ print $1 }'`
#
cat 2dx_merge_dirfile.dat | sort | ${bin_2dx}/2dx_statistics_file.exe > SCRATCH/2dx_statistics_file.dat
#
# This memorizes the current merge directory under the variable "olddir":
set olddir = $PWD
# This translates the list of directories to work on into one single long line:
cat 2dx_merge_dirfile.dat | tr "\n" " " > SCRATCH/2dx_merge_dirfile_oneline.dat
set dirlist = "`cat SCRATCH/2dx_merge_dirfile_oneline.dat`"
#
#############################################################################
${proc_2dx}/linblock "get_max_defocus.py - to extract defocus and TANGL maxima."
#############################################################################
#
${app_python} ${proc_2dx}/get_max_defocus.py 2dx_merge_dirfile.dat SCRATCH/get_max_defocus.out
set tmp_defocus_min = `cat SCRATCH/get_max_defocus.out | awk '{ s = $1 } END { print s }'`
set tmp_defocus_max = `cat SCRATCH/get_max_defocus.out | awk '{ s = $2 } END { print s }'`
set tmp_TANGL_max = `cat SCRATCH/get_max_defocus.out | awk '{ s = $3 } END { print s }'`
#
set tmp_defocus_min = `echo ${tmp_defocus_min} | awk '{ s = $1 / 10000.0 } END { printf ( "%16.3f", s ) }'`
set tmp_defocus_max = `echo ${tmp_defocus_max} | awk '{ s = $1 / 10000.0 } END { printf ( "%16.3f", s ) }'`
set tmp_overall_R_factor = `echo ${overall_R_factor} | awk '{ s = $1 * 100.0 } END { print s }'` 
set tmp_overall_weighted_R_factor = `echo ${overall_weighted_R_factor} | awk '{ s = $1 * 1.0 } END { print s }'`
set tmp_nominal_tilt_angle = `echo ${tmp_TANGL_max} | awk '{ s = int ( $1 / 5 ) * 5 } END { printf ( "%10.0f", s ) }'` 
set realcell_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = $1 } END { printf ( "%10.1f", s ) }'`
set realcell_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = $2 } END { printf ( "%10.1f", s ) }'`
set realcell_gamma = `echo ${realang} | awk '{ s = $1 } END { printf ( "%10.1f", s ) }'`
set realcell_c = `echo ${ALAT} | awk '{ s = $1 } END { printf ( "%10.1f", s ) }'`
#
set tmp_overall_phase_residual_2D = `echo ${overall_phase_residual_2D} | awk '{ s = $1 } END { printf ( "%10.1f", s ) }'`
set tmp_overall_weighted_phase_residual_2D = `echo ${overall_weighted_phase_residual_2D} | awk '{ s = $1 } END { printf ( "%10.1f", s ) }'`
#
#
if ( ${merge_modus} == '3D' ) then
  #
  if ( ${merge_alsoevenodd} == "y" ) then
    # #############################################################################
    # ${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for EVEN images"
    # #############################################################################    
    # echo "# IMAGE: APH/latfitted_even.hkl <HKL: Input HKL EVEN [H,K,L,A,PHI,FOM]>" >> LOGS/${scriptname}.results
    # ${bin_2dx}/2dx_processor.exe --hklin APH/latfitted_even.hkl --mrcout volume_even.map -X ${realcell_a} -Y ${realcell_b} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME}
    echo "# IMAGE: volume_even.map <MAP: final map 3D EVEN images>"  >> LOGS/${scriptname}.results
    #
    echo "<<@progress: +5>>"
    # #############################################################################
    # ${proc_2dx}/linblock "2dx_processor.exe - to transform HKL file into volume for ODD  images"
    # #############################################################################    
    # echo "# IMAGE: APH/latfitted_odd.hkl <HKL: Input HKL ODD  [H,K,L,A,PHI,FOM]>" >> LOGS/${scriptname}.results
    # ${bin_2dx}/2dx_processor.exe --hklin APH/latfitted_odd.hkl  --mrcout volume_odd.map  -X ${realcell_a} -Y ${realcell_b} -Z ${ALAT} --gamma ${realang} --res ${RESMAX} -s ${SYM_NAME}
    echo "# IMAGE: volume_odd.map <MAP: final map 3D ODD images>"  >> LOGS/${scriptname}.results
    echo "<<@progress: +5>>"
    #
    #############################################################################
    ${proc_2dx}/linblock "2dx_correlator.exe - for FSC"
    #############################################################################  
    ${bin_2dx}/2dx_correlator.exe --vol1 volume_even.map --vol2 volume_odd.map -R ${RESMAX} --bins 100 --fsc even_odd_fsc.dat
    #
    ${app_python} ${proc_2dx}/plotFSC.py even_odd_fsc.dat PS/even_odd_fsc.ps
    #
    echo "# IMAGE: even_odd_fsc.dat <TXT: FSC data between even odd datasets>" >> LOGS/${scriptname}.results
    echo "# IMAGE-IMPORTANT: PS/even_odd_fsc.ps <PS: FSC plot between even odd datasets>" >> LOGS/${scriptname}.results
    #
    echo "<<@progress: +10>>"
  endif  
  #
  #############################################################################
  ${proc_2dx}/linblock "2dx_calc_num_phases - to estimate theoretically possible numbers of reflections"
  #############################################################################  
  #
  \rm -f 2dx_calc_num_phases.out
  #
  ${bin_2dx}/2dx_calc_num_phases.exe << eot
${spcgrp}
${realcell},${realang}
${ALAT}
${tmp_TANGL_max}
${RESMAX}
${zstarrange}
eot
  # 
  if ( ! -e 2dx_calc_num_phases.out ) then
    ${proc_2dx}/protest "ERROR in 2dx_calc_num_phases"
  endif
  set num_reflections_theoretic = `cat 2dx_calc_num_phases.out`
  if ( ${SYM} == 'c222' ) then
    echo "::For symmetry c222 reducing the numbers of reflections by 50%,"
    echo "::to compensate for the systematically absent odd lattice lines."
    set num_reflections_theoretic = `echo ${num_reflections_theoretic} | awk '{ s = $1 / 2.0 } END { print s }'`
  endif
  \rm 2dx_calc_num_phases.out
  set completeness1 = `echo ${num_reflections_theoretic} ${num_reflections_FOM1} | awk '{ s = 100.0 * $2 / $1 } END { printf ( "%10.1f", s ) }'`
  set completeness50 = `echo ${num_reflections_theoretic} ${num_reflections_FOM50} | awk '{ s = 100.0 * $2 / $1 } END { printf ( "%10.1f", s ) }'`
  #
  echo "==============================================================================" > Summary_Statistics_3D.txt
  echo "3D Reconstruction Parameters" >> Summary_Statistics_3D.txt
  echo "==============================================================================" >> Summary_Statistics_3D.txt
  echo "Crystal plane group symmetry:# ${SYM}" >> Summary_Statistics_3D.txt
  echo "Crystal unit cell parameters:# a=${realcell_a}A, b=${realcell_b}A" >> Summary_Statistics_3D.txt
  echo "                             #  gamma=${realcell_gamma}deg, c=${realcell_c}A" >> Summary_Statistics_3D.txt
  echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
  echo "Number of images:# ${tmp_num_lines}" >> Summary_Statistics_3D.txt
  echo "Range of defocus [micrometers]:# ${tmp_defocus_min} ... ${tmp_defocus_max}" >> Summary_Statistics_3D.txt
  echo "IQ range used:# 1 ... ${MergeIQMAX}" >> Summary_Statistics_3D.txt
  echo "Tilt range used [deg]:# 0.0 ... ${tmp_TANGL_max}" >> Summary_Statistics_3D.txt
  echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
  echo "In-plane resolution cut-off [Angstroems]:# ${RESMAX}" >> Summary_Statistics_3D.txt
  echo "Vertical resolution cut-off [Angstroems]:# ${zstarrange_real}" >> Summary_Statistics_3D.txt
  echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
  echo "Number of observed reflections:# ${Reflections_Unique}" >> Summary_Statistics_3D.txt
  echo "Number of observed unique reflections:# ${JREFL}" >> Summary_Statistics_3D.txt
  echo "Number of possible unique reflections in asym. unit:# ${num_reflections_theoretic}" >> Summary_Statistics_3D.txt
  if ( ${latline_algo}x == "0x" ) then
  # The following number is without resolution limit:
  # echo "Number of observed unique reflections in asym. unit without resolution limit:# ${num_reflections_fitted}" >> Summary_Statistics_3D.txt
    echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
    echo "Number obs. uniq. refl. in asym. unit with FOM>1%:# ${num_reflections_FOM1}" >> Summary_Statistics_3D.txt
    echo "Completeness to ${tmp_nominal_tilt_angle}deg tilt and ${RESMAX}A (${zstarrange_real}A vertical)," >> Summary_Statistics_3D.txt
    echo "counting reflections with FOM > 1%:# ${completeness1}%" >> Summary_Statistics_3D.txt
    echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
    echo "Number obs. uniq. refl. in asym. unit with FOM>50%:# ${num_reflections_FOM50}" >> Summary_Statistics_3D.txt
    echo "Completeness to ${tmp_nominal_tilt_angle}deg tilt and ${RESMAX}A (${zstarrange_real}A vertical)," >> Summary_Statistics_3D.txt
    echo "counting reflections with FOM > 50%:# ${completeness50}%" >> Summary_Statistics_3D.txt
    echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
    # echo "Overall phase residual [deg]:# ${overall_phase_residual}" >> Summary_Statistics_3D.txt
    echo "Overall weighted phase residual [deg]:# ${overall_weighted_phase_residual}" >> Summary_Statistics_3D.txt
    # echo "Overall R-factor:# ${tmp_overall_R_factor}%" >> Summary_Statistics_3D.txt
    echo "Overall weighted R-factor:# ${tmp_overall_weighted_R_factor}%" >> Summary_Statistics_3D.txt
    echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
  endif
  set lines = `cat SCRATCH/2dx_statistics_file.dat | wc -l`
  set toomany = `echo ${lines} | awk '{ if ( $1 > 20 ) { s = 1 } else { s = 0 } } END { print s }'`
  if ( ${toomany} == "1" ) then
    head -n 20 SCRATCH/2dx_statistics_file.dat >> Summary_Statistics_3D.txt
    echo "::[list truncated] (total directories: ${lines})"  
  else
    cat SCRATCH/2dx_statistics_file.dat >> Summary_Statistics_3D.txt
  endif
  echo "------------------------------------------------------------------------------" >> Summary_Statistics_3D.txt
  #
  #
  echo "::"
  cat Summary_Statistics_3D.txt | sed 's/ /_/g' | sed 's/#/ /g' | awk '{ printf "%-53s%25s\n", $1, $2 }' | sed 's/_/ /g' > tmp
  mv -f tmp Summary_Statistics_3D.txt
  cat Summary_Statistics_3D.txt | sed 's/^/::/'
  echo "::"
  echo "# IMAGE-IMPORTANT: Summary_Statistics_3D.txt <TXT: Summary_Statistics_3D>" >> LOGS/${scriptname}.results
  echo "# IMAGE: LOGS/2dx_tltplotk.txt <LOG: PLTILTK summary>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: PS/2dx_tltplotk.ps <PS: TLTPLOT file>" >> LOGS/${scriptname}.results
  echo "# IMAGE-IMPORTANT: PS/latline.ps <PS: Lattice lines>" >> LOGS/${scriptname}.results
  #
endif
#
echo "# IMAGE-IMPORTANT: PS/2dx_plotresolution.ps <PS: Resolution Plot>" >> LOGS/${scriptname}.results
#
set CSVfile = 2dx_mergeRefine.console.csv
if ( ! -e ${CSVfile} ) then
  set CSVfile = SCRATCH/2dx_mergeRefine.console.csv
endif
if ( -d ../../TILES ) then
  set in_tile = 1
else
  set in_tile = 0
endif
if ( -e ${CSVfile} && ${in_tile} == "1" ) then
  ${app_python} ${proc_2dx}/plotTileBeamTilt.py ${CSVfile} 5 ${CSVoutfile}
  #    
  echo "# IMAGE: ${CSVfile} <TXT: Tile refinement results>" >> LOGS/${scriptname}.results
  echo "# IMAGE: ${CSVoutfile} <PS: Tile refinement plot>" >> LOGS/${scriptname}.results
endif
#
#----------------------------------------------------------------------------------------------------
${proc_2dx}/linblock "2dx_quality_evaluator.exe - to calculate phase residual and fsc of merged data"
#---------------------------------------------------------------------------------------------------
#
set cellx = `echo ${realcell} | cut -d\, -f1`
set celly = `echo ${realcell} | cut -d\, -f2`
#
echo "# IMAGE: APH/latlines.dat <Latline after prescal [H,K,Z,A,P,SigA,SigP,IQ]>" >> LOGS/${scriptname}.results
${bin_2dx}/2dx_quality_evaluator.exe --hkzin APH/latlines.dat --nx ${cellx} --ny ${celly} --nz ${ALAT} --gamma ${realang} --res ${RESMAX} --residuals residuals_merged.txt --fsc fsc_merged.txt  
#
${app_python} ${proc_2dx}/plotPhaseResiduals.py residuals_merged.txt PS/residuals_merged.ps
${app_python} ${proc_2dx}/plotFSC.py fsc_merged.txt PS/fsc_merged.ps
#
echo "# IMAGE-IMPORTANT: residuals_merged.txt <TXT: Phase residuals of merged data>" >> LOGS/${scriptname}.results
echo "# IMAGE-IMPORTANT: PS/residuals_merged.ps <PS: Merged phase residual plot>" >> LOGS/${scriptname}.results
# echo "# IMAGE-IMPORTANT: fsc_merged.txt <TXT: FSC of merged data>" >> LOGS/${scriptname}.results
# echo "# IMAGE-IMPORTANT: PS/fsc_merged.ps <PS: Merged FSC plot>" >> LOGS/${scriptname}.results
#
echo "<<@progress: 100>>"
${proc_2dx}/linblock "Normal End."

exit
# for GUI:
python ${proc_2dx}/get_max_defocus.py 2dx_merge_dirfile.dat SCRATCH/get_max_defocus.out
